/*
 * EpsilonGrid.cpp
 *
 *  Created on: Mar 29, 2012
 *      Author: marchi
 */

#include "EpsilonGrid.h"

namespace EpsilonNS {
	Matrix EpsilonGrid::zero=Matrix().Unit();

	bool EpsilonGrid::operator()(Polarization & P, Field & Ef, Matrix & Rot,  Voronoi * v){
		try{if(!nnx || !nny || !nnz) throw "Grid not initialized for EpsilonGrid ";}
		catch(const char * s){std::cout << s << std::endl;return false;}
		Dvect P_t;
		P_t=P.Integrate();
		Dvect pp,ee;
		for(unsigned int i=0;i<nnx;i++){
			for(unsigned int j=0;j<nny;j++){
				for(unsigned int k=0;k<nnz;k++){
					pp=P[i][j][k];
					ee=Ef[i][j][k];
					D[i][j][k]+=pp % P_t;
					G[i][j][k]+=ee % P_t;
					for(int o=0;o<DIM;o++) M[i][j][k][o]+=P[i][j][k][o];
					for(int o=0;o<DIM;o++) E[i][j][k][o]+=Ef[i][j][k][o];
				}
			}
		}
		M_avg+=P_t;
		TotalCount++;
		return true;
	}


	EpsilonGrid::EpsilonGrid() {

		try{
			if(nnx && nny && nnz) Allocate();
			else throw " Warning Allocation is differred !!!";
		}
		catch(const char * s) {
			std::cout << s << std::endl;
		}
	};
void EpsilonGrid::Rdf(const double dx, const double cut){
		vector<IndexedMatrix> Rdf;
		int nrdf=static_cast<int> ((cut/dx)+1);
		vector<HistData>  rdf(nrdf);
		try{
			if(!is_COset) throw " Must provide Metric ";
			if(!is_xRefset) throw " Must provide reference point ";
		}
		catch(const char * s){
			std::cerr << s << std::endl;
			exit(1);
		}
		Dvect xa;
		const Dvect & x=xRef;
		xa[XX] = oc[XX][XX] * x[XX] + oc[XX][YY] * x[YY] + oc[XX][ZZ] * x[ZZ];
		xa[YY] = oc[YY][XX] * x[XX] + oc[YY][YY] * x[YY] + oc[YY][ZZ] * x[ZZ];
		xa[ZZ] = oc[ZZ][XX] * x[XX] + oc[ZZ][YY] * x[YY] + oc[ZZ][ZZ] * x[ZZ];
		int nx0=(int) nnx;
		int ny0=(int) nny;
		int nz0=(int) nnz;
		double rx,ry,rz;
		double rnx=1.0/ (double)(nx0);
		double rny=1.0/ (double)(ny0);
		double rnz=1.0/ (double)(nz0);
		double dist;
		Dvect xc,xd;
		for(int i=0;i<nx0;i++){
			rx=static_cast<double> ((i)*rnx);
			for(int j=0;j<ny0;j++){
				ry=static_cast<double> ((j)*rny);
				for(int k=0;k<nz0;k++){
					rz=static_cast<double> ((k)*rnz);
					xd[XX]=rx-xa[XX];
					xd[YY]=ry-xa[YY];
					xd[ZZ]=rz-xa[ZZ];
					xd[XX]=xd[XX]-rint(xd[XX]);
					xd[YY]=xd[YY]-rint(xd[YY]);
					xd[ZZ]=xd[ZZ]-rint(xd[ZZ]);

					xc[XX]=co[XX][XX]*xd[XX]+co[XX][YY]*xd[YY]+co[XX][ZZ]*xd[ZZ];
					xc[YY]=co[YY][XX]*xd[XX]+co[YY][YY]*xd[YY]+co[YY][ZZ]*xd[ZZ];
					xc[ZZ]=co[ZZ][XX]*xd[XX]+co[ZZ][YY]*xd[YY]+co[ZZ][ZZ]*xd[ZZ];
					dist=sqrt(xc[XX]*xc[XX]+xc[YY]*xc[YY]+xc[ZZ]*xc[ZZ]);
					if(dist < cut){
						Euler Eul(xc);
						Dvect p=M[i][j][k];
						Dvect e=E[i][j][k];
						Dvect LocP,LocE;
						bool bRot=true;
						if(bRot){
							LocP=Eul.Rotate(p);
							LocE=Eul.Rotate(e);
						}else{
							LocP=p;
							LocE=e;
						}
						int h=static_cast<int> (dist/dx);

						rdf[h].e+=LocE;
						rdf[h].p+=LocP;
						rdf[h].idx++;
					}

				}
			}
		}
		for(unsigned int i=0;i<nrdf;++i){
			IndexedMatrix tmp;
			if(rdf[i].idx){
				tmp.p=rdf[i].p/static_cast<double> (rdf[i].idx);
				tmp.e=rdf[i].e/static_cast<double> (rdf[i].idx);
				tmp.i=i;
				Rdf.push_back(tmp);
			}
		}
		cout << "@  s1" << " legend " << "\"" << "My comm "<< "\"" << std::endl;
		for(unsigned int i=0;i< Rdf.size();i++){
			Dvect LocP=Rdf[i].p;
			Dvect LocE=Rdf[i].e;
			cout << std::setprecision(6) << std::setw(13) << std::fixed << Rdf[i].i*dx*10;
			cout << std::setprecision(6) << std::setw(13) << std::fixed << LocP[XX];
			cout << std::setprecision(6) << std::setw(13) << std::fixed << LocP[YY];
			cout << std::setprecision(6) << std::setw(13) << std::fixed << LocP[ZZ];
			cout << std::setprecision(6) << std::setw(13) << std::fixed << LocE[XX];
			cout << std::setprecision(6) << std::setw(13) << std::fixed << LocE[YY];
			cout << std::setprecision(6) << std::setw(13) << std::fixed << LocE[ZZ];
			cout << std::endl;
		}
		cout << "&\n";
	}

	EpsilonGrid::~EpsilonGrid() {Deallocate();}
	void EpsilonGrid::WriteIt(std::ofstream & fout) const {
		if(!Parallel::comm->Get_Rank()){
			fout.write((char *) &nnx, sizeof nnx);
			fout.write((char *) &nny, sizeof nny);
			fout.write((char *) &nnz, sizeof nnz);
			fout.write((char *) &co, sizeof co);
			fout.write((char *) &co, sizeof oc);
			fout.write((char *) &xRef, sizeof xRef);
		}
		Parallel::comm->Barrier();

		array3<Matrix> D1;
		array3<Matrix> G1;
		array3<Dvect>  M1;
		array3<Dvect>  E1;
		D1.Allocate(nnx,nny,nnz);
		G1.Allocate(nnx,nny,nnz);
		M1.Allocate(nnx,nny,nnz);
		E1.Allocate(nnx,nny,nnz);
		D1=D;
		G1=G;
		M1=M;
		E1=E;
		Dvect M1_avg=M_avg;
		int ntot=TotalCount*Parallel::comm->Get_Size();
		double fact0=1.0/static_cast<double> (ntot);
		M1_avg*=fact0;
		for(unsigned int i=0;i<nnx;i++){
			for(unsigned int j=0;j<nny;j++){
				for(unsigned int k=0;k<nnz;k++){
					D1[i][j][k]*=fact0;
					G1[i][j][k]*=fact0;
					M1[i][j][k]*=fact0;
					E1[i][j][k]*=fact0;
				}
			}
		}
		int Tcount=TotalCount;
		Parallel::comm->ReduceSum(&Tcount,1);
		if(!Parallel::comm->Get_Rank())	fout.write((char *) &Tcount, sizeof Tcount);
		Parallel::comm->ReduceSum(&D1[0][0][0][0][0],DIM*DIM*nnx*nny*nnz);
		Parallel::comm->ReduceSum(&G1[0][0][0][0][0],DIM*DIM*nnx*nny*nnz);
		Parallel::comm->ReduceSum(&M1[0][0][0][0],DIM*nnx*nny*nnz);
		Parallel::comm->ReduceSum(&E1[0][0][0][0],DIM*nnx*nny*nnz);

		if(!Parallel::comm->Get_Rank())	{
			fout.write((char *) &D1[0][0][0][0][0], (sizeof D1[0][0][0][0][0])*DIM*DIM*nnx*nny*nnz);
			fout.write((char *) &G1[0][0][0][0][0], (sizeof G1[0][0][0][0][0])*DIM*DIM*nnx*nny*nnz);
			fout.write((char *) &M1[0][0][0][0], (sizeof M1[0][0][0][0])*DIM*nnx*nny*nnz);
			fout.write((char *) &E1[0][0][0][0], (sizeof E1[0][0][0][0])*DIM*nnx*nny*nnz);
		}
		Parallel::comm->Barrier();
	}
	void EpsilonGrid::ReadIt(std::ifstream & fin){
		Matrix co0,oc0;
		fin.read((char *) &nnx, sizeof nnx);
		fin.read((char *) &nny, sizeof nny);
		fin.read((char *) &nnz, sizeof nnz);
		fin.read((char *) &co0, sizeof co0);
		fin.read((char *) &oc0, sizeof oc0);
		fin.read((char *) &xRef, sizeof xRef);
		fin.read((char *) &TotalCount, sizeof TotalCount);
		Deallocate();
		is_xRefset=true;
		is_COset=true;
		setMetric(co0);
		Allocate();
		fin.read((char *) &D[0][0][0][0][0], (sizeof D[0][0][0][0][0])*DIM*DIM*nnx*nny*nnz);
		fin.read((char *) &G[0][0][0][0][0], (sizeof G[0][0][0][0][0])*DIM*DIM*nnx*nny*nnz);
		fin.read((char *) &M[0][0][0][0], (sizeof M[0][0][0][0])*DIM*nnx*nny*nnz);
		fin.read((char *) &E[0][0][0][0], (sizeof E[0][0][0][0])*DIM*nnx*nny*nnz);
	}


} /* namespace EpsilonNS */
